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HALF-SPACE KINETIC EQUATIONS WITH GENERAL BOUNDARY CONDITIONS 


QIN LI, JIANFENG LU, AND WEIRAN SUN 


Abstract. We study half-space linear kinetic equations with general boundary conditions that consist of 
both given incoming data and various type of reflections, extending our previous work [LLS14] on half¬ 
space equations with incoming boundary conditions. As in [LLS14] , the main technique is a damping 
adding-removing procedure. We establish the well-posedness of linear (or linearized) half-space equations 
with general boundary conditions and quasi-optimality of the numerical scheme. The numerical method 
is validated by examples including a two-species transport equation, a multi-frequency transport equation, 
and the linearized BGK equation in 2D velocity space. 


1. Introduction 


In this paper we propose an efficient numerical method for linear half-space kinetic equations with general 
boundary conditions 

fida:/+ Cf = 0 , in(0,oo)xV, 

/h>o = ^ uh<o) ’ on X = 0 , 

where the density function f{x,v) G K"* with m > I for x G [0, c») and v = = {p.,V 2 ,- ■■ ,Vd) G V. 

Typical examples for the velocity space V include the whole space as in the case of the Boltzmann 
equation, and V = [—l,l] ss in the case of the transport equation. By allowing higher-dimensions in / 
and V, we include multi-species models and models with multi-dimensional velocity variables such as the 
linearized Boltzmann and linearized BGK equations. The setup also includes the multi-frequency case where 
the frequency variable can be treated as an index for multi-species after discretization. 


The operator £ in (l.l) is a linear operator, examples of which include the scattering operator in the 


linear transport equations, the collision operator in the linearized Boltzmann equations and the linearized 
BGK equation. The operator JC is the boundary operator which characterizes various types of reflections 
at the boundary. Two classical examples for the reflections are the diffuse and specular reflections. It will 
be discussed in details in Section |2.2| that our method applies to a general class of boundary operators 
including Maxwell boundary condition (linear combination of the diffuse and specular reflection), bounce- 
back reflection, and also the more general (linearized) Gercignani-Lampis boundary condition. 

It is well known that to ensure the well-posedness of equation HD), one needs to prescribe suitable 


boundary conditions at a; = oo. The precise conditions were hrst formulated in GGS 88 for the linearized 


Boltzmann equations with prescribed incoming data. This type of well-posedness result has been extended to 
general linear/linearized half-space equations and weakly nonlinear half-space equations with both incoming 
and Maxwell boundary conditions (see e.g., G 0 IO 8 STll UYY 03| ) and also to discrete Boltzmann equation 
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with general boundary conditions Ber08 BerlO . This is also the setting that we use for developing numerical 


methods. Now we briefly explain the details of the formulation of the boundary condition at infinity. Denote 
the null space of £ as Null£ which is assumed to be finite-dimensional. Let V be the L^-projection operator 
onto Null£ and as the corresponding orthogonal projection operator such that 

IP : L^(dCT) ^Nulir, V^=l-V. 

Define the operator Vi : Null£ —> Null£ as 

Vj=V{fif) for/GNull£. 

It is clear that Vi is a symmetric operator on a finite-dimensional space, and hence all its eigenvalues are 
real. Denote the eigenspaces of Vi associated with positive, negative, and zero eigenvalues as 
respectively. Then Null£ is decomposed as 

Null/: = iL+0 iJ-© . 


Using these notations, we prescribe the boundary conditions at a; = oo in a similar way as in CGS88 such 
that 

lim f G H+ . 

The complete form of the kinetic equation considered in this paper reads 

flda:f + Cf = 0, 




I /i>0 
lim f G 

X—^OO 


( 1 . 2 ) 




More specific assumptions regarding C and K, to guarantee the well-posedness of (1.2) will be discussed in 
Section [2l 

Half-space equations with general boundary conditions are frequently encountered in electric propulsion for 


satellites GKOS and photon transport in solid state devices HM14 HRBIO , among many other applications. 


The standard treatment of this type of equations is the Monte Carlo method HRBIO . There are also special 


cases where analytical solutions are possible HM14 


In LLS14 we developed a direct systematic method to solve half-space equations in the case of pure 


incoming boundary condition (when /C = 0). There are also other direct numerical approaches for this 
case proposed in Cor90 GK95| . Compared with our approach, these methods suffer from severe Gibbs 
phenomena and lack of error analysis or systematic strategy to reduce numerical errors. We also note that 


the method for linearized discrete equations in Ber08 can be applied to solve the continuous half-space 


equation by approximating it using discrete velocity models. Unlike Ber08 which focuses on the analysis 


of the discrete model, our goal here is to approximate the solutions to the continuous half-space equation 
using a spectral type method with convergence analysis. 

The present work extends our previous method to the case when various reflections are involved. The main 
difficulties that we need to overcome are the degeneracy of £, the derivation of a proper weak formulation 
involving ]C, and the fact that the boundary conditions at a; = oo are part of the solution instead of being 
prescribed. To this end, we apply similar procedure proposed in LLS14 , which combines and extends 


the ideas of even-odd decomposition ES12 and a damping adding-removing procedure UYY03 Gol08 


More specifically, we first modify £ by adding damping terms to it (see Section 2.1). This will remove the 


degeneracy of £ and ensure that the end-state of the damped solution at a; = oo is zero. Both analysis and 
numerical schemes are then performed on the weak formulation of the damped equation, which is derived 


by applying an even-odd decomposition with mixed regularity ES12 of the (damped) solution /. One 
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important advantage of the even-odd decomposition is that it leads to a natural way of constructing a 
family of basis functions that captures the possible jump discontinuity of the solution at /r = 0, by the odd 
extension of the basis functions constructed for positive ii. This discretization of the velocity space based 
on even-odd decomposition turns out to be equivalent to the double-method developed in the literature 
of solving neutron transport equations, see e.g., Tho63 . We also comment that the appearance of the 


boundary operator K, introduces extra difficulty into formulating the weak form of the half-space equation. 
The difficulty comes from the fact that only the even part of the solution /+ has enough regularity to define 
a trace on the boundary. Our main idea here is to use the properties assumed for K, in Section 2 to represent 
the odd part /“ on the boundary in terms of /“*■. 

Our numerical method is spectral in nature: we apply Galerkin approximations to the weak formulation 
and use Babuska-Aziz lemma to show that the damped equation is well-posed and the finite-dimensional 
approximation is quasi-optimal. It will be clear that the damping plays a crucial role here. Finally, we 
make use of the linearity and use proper superposition of certain special solutions to the damped equation 
to recover the original undamped solution. 

A by-product of the above procedure is that we obtain a unified proof for the well-posedness of the half¬ 
space equations with general boundary conditions. This well-posedness theory is general enough to include 
multi-species and multi-dimensional (in velocity) half-space equations. 

The layout of the paper is as follows. In Section 2 we explain all the assumptions for the linear operator C 
and the boundary operator JC. In Section 3 we prove the well-posedness of the half-space equation using the 
damping adding-removing procedure. In Section 4 we show three numerical examples which cover the three 
cases of multi-species, multi-frequency transport equations and a multi-dimensional (in velocity) linearized 
BGK equation. 


2. Main Assumptions for C and JC 

In this section we collect the conditions on the linear operator C and the boundary operator 1C. 
Notation. In this paper we denote 

(/, 5 ) = [ f-gda, and (f, g) = f [ f-gdadx, 

' ' « Jv ' ' Jr Jw 

where du is a measure in the velocity space. Throughout this paper we assume that the measure dcr is 
symmetric with respect to /r. 

2.1. Main Assumptions for C. In this subsection we state the general assumptions for the collision 
operator C. First, define the weight function (attenuation coefficient) 

a(u) = (I + (2.1) 

for some 0 < kq < 1- The first four basic assumptions for the linear operator C are as follows: 

(PLl) C : 'D{C) —)• (L^(dcr))™ is self-adjoint with its domain 'D{C) given by 

V{C) = {/e {L\da)r\ a{v)fe {L\da)r} C {L\da)r , 

where a{v) is defined in ( |2.1[ ). Such space arises naturally for linear/linearized collision operator 
since in many cases C has the structure as 

C = a{v)I + Cl , 

where Ci is a bounded or even compact operator. 
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(PL2) C : (L^(ad(T))™ —?► dc))’” is bounded, that is, there exists a constant Co > 0 such that 

(L^(a dtr))’^ 

(PL3) Null£ is finite dimensional and Null£ C (£P(dCT))"* for all p € [l,oo). 

(PL4) £ is nonnegative: for any / G (£^(ado-))™. 


Cf 




<Co 


/ • £/dcr > 0. 


( 2 . 2 ) 


Assumptions (PL1)-(PL4) are general enough to include many classical models such as the linearized 
Boltzmann operators (around Maxwellians) with hard-potentials, the linearized BGK operator, and linear 
transport operators for single- or multi-species. In fact, these classical operators satisfy an even stronger 
coerciveness property: 


/ • £/ dcr > Co 


p-L/ 


(L^ (a d(7))^ 


(2.3) 


where recall that V'^ = I — P and V : (L^(d(T))'" —>■ Null£ is the projection onto Null£. 

We need one last essential assumptions on the coercivity of a damped version of £ on the whole (£^( dcr))™ 
but not just (Null £)■*■. To properly explain this assumption, we introduce several definitions related to the 
null space of £. Recall that Vi : Null£ —>■ Null£ is the operator given by 


'Piif) = ^(£/) for any / G Null£ . 


Note that Vi is a symmetric operator on the finite dimension space Null£. Therefore, its eigenfunctions 
form a complete basis of Null£. Denote ,H^ as the eigenspaces of Vi corresponding to positive, 

negative, and zero eigenvalues respectively and denote their dimensions as 

= dim H '^, i/_ = dim H ~, i/q = dim . 


Let Al_|__i, Al_j,Xo,fe be the associated orthornormal eigenfunctions with 1 < i < u_|_, 1 < j < iz_, and 
1 < fc < uq. Note that if any of u±, vq is equal to zero, we simply do not have any eigenfunction associated 
with the corresponding eigenspace. By definition, these eigenfunctions satisfy 

(^r, 7 ) 5 — 0 if T ^ T Or y ^ y , 

Xo^k)y = o, x+^i)^>o, (/iAr_j, Ar_j)^<o, 

where r G {-I-,— ,0}, y G {i,j,k}, 1 < i < u_|_, 1 < j < u_, and 1 < k < i/q. 

Our method relies on full coercivity of the collision/scattering operator on (£^(adcr))™ instead of the 
partial one in (2.3) on (Null£)^. Hence, instead of working directly with £, we add in the damping terms 
on the modes in Null£ and define the damped linear operator Cd as 

1 '+ 


Cdf=Cf-\- a / 






(2.4) 


/i£ ^{pXo,k) ^{pXo^k), , 

k=l ^ k=l “ 

where a > 0 is some constant damping coefficient to be determined later. The motivation of defining Cd 
in such a form is as follows: the operator £ normally will provide bounds for the orthogonal component of 
/ in (Null £)'*■. With the added damping terms to dissipate the modes in Null£, we expect that Cd will 
satisfy certain full coercivity condition on (£^(d(T))™. On the other hand, this added damping effect can 
be eventually removed using linearity of the equations. The precise assumption of £ regarding its coercivity 
states 
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(PL5) There exist two constants a, fJo > 0 such that the damped operator Cd satisfies 

2 

(L^(a do'))’^ 

for any /G {L^{aAa))”^. 


f ■ Cdf dcr > (To 


/ 


(2.5) 


It will be shown in Lemma |4. 2 [ that the coercivity condition in (2.31 combined with the form of Cd in (2.4) 
implies (PL5), and hence (PL5) is a natural assumption for many examples. 

2.2. Main Assumptions for 1C. In this part we specify conditions for the boundary operator JC. These 
conditions are stated in rather general forms and are satisfied by a large class of boundary operators. Recall 
that we have denoted v = (/i,u) = (/r, U 2 , • • • ,Vd)- Denote the incoming and outgoing parts of the velocity 
space as 

V+ = {u = (/X, U) I /r > 0} , and V_ = {v = \ fj. < 0} . 

We consider the general case where the boundary operator /C consists of various types of reflections in the 
sense that there exists a coefficient G [0,1) and a scattering kernel kr (which is a positive measure) such 
that 


JC = arJCr , [JCrf]{v) = / kr{v, v') f {v') da{v') for u G V+. 




( 2 . 6 ) 


The main assumption for such 1C is 
(PK) The reflection operator JCr satisfies that 


A* 


' /i>0 


JCrf 


dcr < 


ImI 


' /i<0 


/ 


dcr . 


(2.7) 


There is a large family of reflection boundary operators ICr that satisfy (PK). In the literature, the 
reflection boundary operator for nonlinear kinetic equations of a single species is usually written as 


[jerF]{v) = - [ v')F{v') du', ve V* 

Jd'<0 


for some scattering kernel R. If we consider the linearization around the equilibrium state such that 

F = M + Mf, 

then the linearized version has the form 

M-i(u) 


[JCrf]{v) = 




> /i'<0 


fi'R{v, v')f{v')M{v') du' 


( 2 . 8 ) 


We show in the following lemma that as long as JCr satisfies the classical normalization and reciprocity 
conditions, then the main assumption (PK) holds: 

Lemma 2.1. Suppose M is a scalar equilibrium state and dcr = M dv where du is the Lebesgue measure. 
Suppose JCr has the form as in ( |2.8[ ). If R satisfies the normalization and reciprocity conditions: 

(2.9) 

( 2 . 10 ) 


\qi'\M(v')R{v, v') = \pl\M{v)R{—v', —v ), v G V"*" , v' G 
j R{v, v') du = 1, v' C V“ , 

J u>0 


then (PK ) holds. 
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Proof. Note that an immediate variation of (2.10) is 


f R{—v',—v)dv' = l, vGY~'~. 

J fi'<0 


( 2 . 11 ) 


By the definition of we have 


\ICrf\ M{v) dw = 




'/J>0 
< 


' ti>0 

\fi'\R{v, v')p{v')M{v') di; 


f R{v,v')f{v')M{v') dv'^ di; 

Jn'<0 / 

■ w\ 




/i>0 \Jfj,'<0 


'<0 ^ 


-R{v,v')M{v')M ^(z;)dz;M di; 


' fi>0 \J fj, 


\^j,'\R{v,v')f^{v')M{v')dv'j f f R{—v',—v)dv'j dw 

ii'<0 / \JfL'<0 J 


|| 2 . 11 [ 


//x>0 Jfi'<0 


\lx'\R{v, v')f^{R)M{R) diR d^; 


(|/r'|/^(i;')M(v')) f f R{v,v') dv\ di;' f \n'\f'^{v')M{v') dv'. 

u '<0 


The condition (PK) follows as da = M dri in this case. 


□ 


Examples that satisfy (2.9) and (2.10) include 

• the specular reflection condition where i?(u, v') = 5(^ + ^i')5{v — v')\ 

• the bounce-back condition where R(y, v') = S{v + v')] 

• the pure diffuse condition for BGK or linearized Boltzmann equation where R(v, v') = 

• convex combinations of the above three; and more generally, 

• the (linearized) Cercignani-Lampis collision operator with R given by 


(27r) 2 


_ Ml 

-e 2 


R{v,v') = 


1 


■ exp 


-I- (1 - Qn)(A^')' 

2Qr7, 


where 0 < an < 1, 0 < < 2, and 


exp 


cos (p 


' - (1 - at)v' 

2at{2 - at) 


Jo 


( a /1 - anjilJ 
\ 


d^ . 


Hence our method applies to all of these classical cases for single species. 

Remark 2.1. In all of our numerical examples in Section]^ we use either the Dirichlet boundary condition 
with given incoming data or the classical Maxwell boundary condition where 

JC = adlCd + asJCs 

with the accommodation coefficients ad, as satisfying 

ad,as>0, 0<ad-l-as<l. 

The two operators ICd,ICs are the diffuse and specular reflection operators respectively. In terms of the 
notation in (2.6), we can choose in this case 

_ I ir — ir I ir 

O^r — “r CXg 5 A^r — Rd Rs • 

CX’P Ctf 

Since /Cg automatically satisfies (PK) with an equal sign, we only need to check in each numerical example 
that ICd satisfies (PK) as well. 

Below we state two essential consequences of assumption (PK), which will guarantee the well-posedness 
of the half-space equation and provide the foundation for the numerical scheme. 
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Lemma 2.2. Suppose K, satisfies fP'K.). Then the half-space equation (1.2) has at most one solution f € 
C([0,oo);(L2(|^|da)n. 

Proof. By the linearity of the equation we only need to prove that if /i = 0 in the boundary condition of ( |1.2[ ) , 
then the only solution to (1.2) is zero. By the non-negativity of C, we have /y/i|/ | {x, •) dcr is decreasing 
in X. Since fao G H~^ 0 H^, we have /y/r|_^|^dCT > 0. Hence /y ^ |/]^(x, •) da > 0 for all a; > 0. In 
particular this shows 

[ ^i\ff{0,■)da>0. (2.12) 

Jv 

By assumption (PK), at the boundary x = 0 we have 


' /i>0 


dcr = 


' Ai>0 




dcr < 


ImI 


> /i<0 


dcr. 


Therefore, 




/ 


(0, •) dcr < - (1 - a^) / 1^1/ (0,-)dcr<0. 

J ^<0 


By \2.12\ and that 1 — > 0, we deduce that 

f fi f (0, •) dcr = 0 , and hence 

J fi<0 


> fji>0 


(0,-)dcr = 0. 


Therefore, at a; = 0 we have /(O, •) = 0. By the uniqueness of solutions to (1.2) with only the incoming data 
CGS88| (that is, ar = 0), we have that (1.2) has at most one solution. 


□ 


Remark 2.2. The assumption that < 1 in (2.6) is necessary for the uniqueness of the solution. For 
example, if Or = 1 or = 1 in the boundary operator K, for the linear transport equation in (4.9) 

considered in our numerical examples, then any multiple of Xq is a solution to the half-space equation with 
zero incoming data. 

The second consequence of assumption (PK) is 

Lemma 2.3. Suppose the measure dcr in the velocity space is symmetric with respect to p,. Define the 
operator 1C : (L^(/rly>o da))™ —?► (L^(/ily>o dcr))™ such that 

X = aXr , (2.13) 

where Xr is defined as 

{Xrf)ip,v)= f kr{ip,v'),{-p',v'))f{p',v)da{p',v'), Ai>0, 

J /x'>0 

where k^. is the reflection kernel of Xr ■ Note that we have reflected the p' component of v' in the kernel. 
Then 

(a) l-\-X is invertible on (L^(p,ly>o dcr))™. 

(b) There exists a constant fli > 0 such that the operator (I + /C)“^(I — K.) satisfies that 

(pfl (I + (^)-'(I - X)f) > /?1 (/i/, /) (2.14) 

for any / G (L2(^ly>o dcr))™. 


M>0 


Proof, (a) Denote gi{p',v') = f{—p',v'). Then ICrf = Xrfii by the symmetry of dcr with respect to p. 
Hence, 




' A^>0 


Xrf 


dcr = 


p \Xrgi I dcr < 


|gi I dcr = / p 


' M>0 


' A^<0 


> fL>0 


dcr. 
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Therefore, 




If ll£((L2(;,l^^od<T))™) < “r < 1- 

This shows I + /C is invertible on (L^(//l^>o dcr))"^. Furthermore, we have the bound 

_ _ ||2.15 

||(I + K.) ^ ^ (l “ IF IL((L2(^l^>od^))™)) 

(c) Denote 52 = (1 + ^)”^/G (L^(/rl^>o do-))™. Then 

(m/, (l + ]C)-\l-lC)f)^^^ = (m(I + -^)52, (I--^)52)^>o = (W2, 52)^>o- IC 92 ) 

2 II_II 2 

IIS 2 \\(L^(^l^^oda-))'^ 11^52 dcr)) 

Ils^ ■ 

Observe that 

/ 

(L 2 (^ 1 ^>q do-))" 

We conclude by combining the previous two estimates such that 

(fif, {l + JC)-\l-IC)f) > {l-al){l + ary 

\ / fi >0 


(2.15) 


Ai >0 


< l|(I + ^) ll£((L2(^i^^„da))-") IF 2 ll(L=(MO>odO)™ - 1152 II (L2(^i^^„ d<T)) 


(L 2 (^ 1 ^>P do-))" 


Hence (2.141 holds with /3i = (l — a^)(l + a^) 

3. Well-posedness 


□ 


In this section we establish the well-posedness of equation (1.2 1 based on the assumptions for C and K, 
in the previous section. The framework is similar to LLS14 : first we add damping terms to C and show 
that the damped equation has a unique solution. This will be achieved by using the Babuska-Aziz lemma. 
Then we show how to recover the solution to the original kinetic equation using suitable superpositions with 
special solutions. 

The damped kinetic equation has the form 

Hda:f+ Cdf = 0 , 


/b>0 = 5 + WbJ. f‘>». 

/ —>■ 0 , as a; —>■ oo, 


(3.1) 


where the damped operator Cd is defined in (2.4). 


3.1. Weak Formulation. In order to show the well-posedness of (3.1), we consider the weak formulation 
of the equation using the even-odd decomposition. Recall that we have denoted v = {v 2 , ■ ■ ■ ,Vd)- For any 
scalar function 9 ( 9 , v), let 9 ^,g~ be its even and odd parts (with respect to 9 ) respectively such that 


9+{9,v) = 


9{n,v) + 9{-ti,v) 


9 {9',v) = 


9{n,v) - 9{-9,v) 


Therefore we have 

5 + ( 5 , v) + g~ { 9 ., v) = 9 ( 9 , v), 
For the vector-valued function /, denote 

T = (/i+, ft, ■■■, ftf, 


9 { 9 , v)-9 i9,v) = g{-9,v). 


r = (/r, fy,---, fyr ■ 
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The solution space for (1.2) and (3.1) is 


r = {/ e iL^idadx)r I e iL^{daidx)r} . 

for some dcri such that the term Jy / • is well-defined. The norm in T is defined as 


/ 


(L2( dcr 


fJ'dxf'' 


{L^{ d(7i dx))” 


(3.2) 


One example of dcr, dui is that dtr = a{v) dn and dcri = ^ dv where a(v) is the attenuation coefficient 
defined in (2.1). For the operator defined in (|4.21 ), we have dui = dcr. 


This type of solution space T with mixed regularity is introduced in 


ES12 


. For a general function / G F, 


the trace of /i/'*' at x = 0 is well-defined while the trace of /i/“ may not. Due to this lack of regularity for 
f~, when deriving the weak formulation we will represent f~ in terms of /+ on the boundary. Recall that 
the boundary condition is given by 


f] 


fi>0 


= h + JC 




Using the even-odd decomposition, we have 


(/^+ r) L>o = h++ ^(r L<o) 


— h Ck.j- 

= h ar 


' A ^'<0 


v), (/r', dcr -h a. 




kr{in,v),{fj.',v'))f {fi')da 










= h + lCf+- jcr 


where K, is defined in (2.13). Note that in order to get the third line we have used that dcr is symmetric 
with respect to fi. Hence, the boundary condition has been reformulated as 

{i + ic)r = h-{i-ic)f+, fi>o. 


For any operator K, that satisfies assumption (PK), we have shown in Lemma 2.3 that I -I- Ai! is invertible on 
(L^(/rly>o da))"*. Thus, f~ is related to /+ as 


r 


I /x>0 


Im >0 ' 


(3.3) 


M >0 




+ (</>) kldf 


Hence when deriving the weak formulation of the half-space, the boundary term at x = 0 becomes 
(m/, = 2 (/r^+, /- = 2 (I + X)-^h )- 2 (m^+, (I + ^)-i(I - K)f+ ) 

Define the bilinear form 

B{fj) = -{r,iJidx$+ 

\ I Xd 

+ 2 (^^+, (i+:^)-i(i-^)/+ 

and let I be the linear functional on (L^(/il^>o dcr))™ such that 

l{$)=2(n$+, {l + ]C)-^h) 

\ / Ai >0 


/i>0 


x—0 


The previous calculations then show that the weak formulation of equation (3.1) has the form 

Bif, $) = l{$) for any ^ G F . 


(3.4) 


(3.5) 


(3.6) 


The main tool that we use to show well-posedness and quasi-optimality is the Babuska-Aziz lemma which 
we recall below: 
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Theorem 3.1 (Babuska-Aziz). Suppose T is a Hilbert space and : F x F 
Let I : T ^ M. be a bounded linear functional on F. 

(a) If B satisfies the boundedness and inf-sup conditions on F such that 

• there exists a constant cq > 0 such that \B{f,g)\ < co||/||r|l5l!r for all f,g G F; 

• there exists a constant i5o > 0 such that 

sup B{f,fi)>6o 


is a bilinear operator on F. 


r ) 


sup B{f,tlj)>6o\ 

bllr = l 


for any tp G T , 
for any f GT 


(3.7) 


for some constant (5o > 0, 
then there exists a unique f GT which satisfies 

Bif,ip) = forany-ipGT. 

(b) Suppose T is a finite-dimensional subspace ofT. If in addition BiT^xYn^R. satisfies the inf-sup 
condition ouTn, then there exists a unique solution fjs; such that 

B{fN,'4’N) = K-fiN), for any ipN G TN . 

Moreover, /tv gives a quasi-optimal approximation to the solution f in (a), that is, there exists a constant 
Ki such that 

||/-/A||r<Ki inf ||/-w||r. 
wGFn 


Now we verify that B{-, •) and l{-) defined in (3.4) and (3.5) satisfy the conditions in Theorem 3.1 


Proposition 3.2. Suppose the measure da in the velocity space is symmetric with respect to fi. Suppose the 
linear operators C satisfies the assumptions (P'L1)-(P'L5) and the boundary operator K. satisfies assumption 
(PK). Then 

(a) the bilinear form S : F x F —> K satisfies the boundedness and inf-sup conditions and the linear 
functional I is bounded on F. Therefore, equation (3.6) has a unique solution f GT. 

(b) Moreover, ^dxf G {if da dx))'^. Thus f is a strong solution to the damped half-space equa¬ 
tion (3.1). 

Proof. For each f G T the proof is done by finding an appropriate test function cpf G T such that B{4>, f) 
satisfies the inf-sup condition: 

2 


BififJ) > Co 

The particular choice oi fif is the same as in 

$1 = f, 


f 


(j)f 


< Cl 


f 


LLS14 


such that (pf = Sifii -\-(j )2 with i5i > 0 large enough and 
^ -l^dxf^ ■ 


(1 -k 

Using such (pf together with the coercivity of the damped operator Cd in (PL5), we have identical estimates 
for the interior terms in B{(pf, f) as in the proof of LLS14, Proposition 3. 1]. M oreover, the positivity of the 

Hence by the same argument 


boundary term ( ptf'^, (I -I- /C) ^(I — A)/'*' ) is guaranteed by Lemma 

u>o 

as m ’ .. . 


2.3 


LLS14 , we have that B satisfies the inf-sup condition. Boundedness of B and I can be shown by direct 


applications of the Cauchy-Schwarz inequality. Thus the weak formulation (3.6) has a unique solution. This 
also implies that the half-space equation (1.2) has a unique solution in the distributional sense. In addition, 
the half-space equation itself shows fid^f = —Tf G (L^(^ dadx))"^ where a is the attenuation coefficient 
defined in (2.1). Hence the full trace of / in L^(|/i| dcr) is well-defined. □ 
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orthonor- 


As in LLS14| we will solve the damped equation ( |3.1[ ) by a Galerkin method. 

Proposition 3.3 (Approximations in K'^). Suppose (v2)---'<p^n}{vd)} 

mal basis of (da) such that 

• V'2n-i(M) V’2n iP-) evcTi in fj, foT any n > 1; 

• (m) e span{^/>J^\ • • • , for each n>l. 

Define the closed subspace Tcls 

r m 2N-\-l K 'I 

^NK = \ g{x,v) G r g{x,v) = {vd) e, \ , 


i—1 k — 1 712,,71d — ^ 


where = (0, • • • , 0,1, 0, • • • , 0)^ is the standard i*^ basis vector with 1 < i < m and G 


Then 

(a) there exists a unique f^K G F such that 

m 2Af+l K 


fNK{x,v)=J 2 J 2 


a 


(b 

k,n 2 ,--- ,nd 


ix)'ipk \T)'ip^^Hv2) ■ ■ ■ i^if]{vd) e*, 


i—1 k — 1 712,■■■ ,7^(1 —1 


which satisfies 


B{fNK,g) = l{g) for every g GTnk 


where B and I are defined in (3.4) and (3.5) respectively. The coefficients {a^ 
that 

Ji) 


k,n 2 ,--- ,nd 


(3.8) 

(3.9) 
(x)} satisfy 


“fc.ns," .nd(’) S C'^[0,oo) n iF^(0, oo), 1 < k < 2N + 1, I < n 2 , - ■ ■ ,n-d < K , 1 <i <m. 

(b) The approximation is quasi-optimal, that is, there exists a constant K 2 > 0 such that 

\\f-fNK\\r<K2 inf \\f-w\\T. 

UieTNK 

Proof. Part (a) and (b) follow directly from the Babuska-Aziz lemma as long as we verify that B{-, •) satisfies 
the inf-sup condition over Fjvic. The only modification is in the choice oi (fj where (j )2 is projected onto Fjvk. 
The proof again follows along the same line to the proof of Proposition 3.2 in LLS14 using the positivity 
of the boundary term guaranteed by Lemma 2.3 □ 


The following Proposition reformulates (3.9) into an ODE with explicit boundary conditions. 

Proposition 3.4. Let 


A = 




(2Af-|-l)x(2Af-|-l) 


Define the 2(d-l-1)-tensors 21 and 25 as 

2l = A(g)/(g)---(g)/(g)J = {AikSn2l2 ■ ■ • ^ndld^pq) (^2N-i-iPxK^X---xK^xm‘^ 


'^kn2"-7idP 


^k^ {t)'^u2 (^'a) • • • {vd) ep, Cd (A)V'lf M ■ ■ ■ -iplf M e,) 


(3.10) 


for 1 < i^k < 2N + 1 < n2, * * * ^nd ^ K, 1 < ^2? * ‘ ‘ fid ^ K, and \ < p^q < m. Then the variational form 

(3.9) is equivalent to the following ODE for the coefficients o!'kn 2 ---nd^^)‘ 


K 


771 2N-\-l 

EE E 

p—l k—1 n2,--- ,nd — l 


Qiih-ldq On-- , 

'^kn2---7l(iP ^^^k7l2---7ld 


Ap) 


K 


m 2N-\-l 

) = EE E 

P—I k—1 7l2,--- ,71d — l 


mih-ldq (P) / -> 

kn2-ndP ^k7i2---7id 


(3.11) 
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together with the boundary conditions at x = 0 : 

N+l 


fc=i 


= 2 


' M >0 


^i{l + lC)-^h■ if. 


2 i,n2,"- ,nd 


U>0 

dfj 


(3.12) 


fori = !,■■■ , iV, 712, • • • , = 1,2, 

feasis function ifff^ ■■■ nw 


,K, and q = I,- - ■ ,d. Here f^K is defined in (3.8 1 at x = 0 and the 




■(«) 

2 i,n 2 ,-" ,n-d 


= M--- i’ffi (Vd) ' 


Proof. Equation (3.11 1 is obtained by choosing the test function g in ( |3.9[ ) as 

g = go{x) ■ ■ ■if^^}{vd)ep 

for each basis function ipf'\tJ-)il^nf {V 2 )--- ll’n] (vd) Gp and for any arbitrary gof) G G( 0 j o®)- The boundary 
condition (3.12) is derived by choosing the test functions as 

g = 9i(x)ip2k M---iplf](vd) Gp, 

for each basis function ip^f ig)ipnf (^ 2 ) • • • V'if (vd) Gp in r NK and for any arbitrary gif) G C).[ 0 , 00 ). □ 

Since the tensors 21, *8 are the same as in [LLS14 , we have that there are mNK‘^~^ positive, mNK‘^~^ 
negative, and generalized eigenvalues of (21,25). Note that there are m{2N + unknowns in 

the ODE system (3.11) and mNK‘^~^ boundary conditions. This is again the correct number of boundary 
conditions for (3.11) to have a unique decaying solution. 

3.2. Recovery. In this part we show the procedures to recover the solution to the original kinetic equa¬ 
tion (1.2). To this end, let / be the solution to the damped equation (3.1 1 . For all 1 < i < hq and 1 < j < 7 +, 
let ffo.zjffoj' be the solution to (3.1) with h = Xo,i ~ and h = ^+7 — /C(X_|_j-j^^^) respectively. 

More explicitly, for each 1 < i < i^o, 


h-dxgo.i + Cdgo,t = 0 , 


9oAp>o- {^o,i - ICiXoA^^A) +^(5o,i|^<o) 


go,i “t 0 , 

and for each 1 < j < 

g^xg+j ~G P-dg+.j — 0 ^ 


g>0, 

clS X —y 


(3.13) 


3 )) (ff+j |p<o) > 


^ > 0 , 
as a; —; 


(3.14) 


9+Ap.>o-\^+d 

9+,j b) 

The key idea is that the damping terms in Cd vanish for a proper linear combination of /, go^s, and goj^s. 
The recovering procedures rely on the uniqueness of solutions to the original kinetic equation (1.2). 

Proposition 3.5. There exists a unique sequence of constants co,i,c+^fc G K for 1 < i < 70 and 1 < j < 7 + 
such that if we define 

70 7 + 

g 'y ^ ^Q igQ i -|- 'y ) Cj^ jgj^ j , 


i=l 


i=l 


(3.15) 
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then 


[t^Xo,ioAf - 9){x)J^=0, (/-5)(a;)^^ =0, (^^lC ^ {hXq^A, {f - g)ix) j ^ = 0 (3.16) 

for all X > 0 and all 1 < < ^q, 1 < i± < 7 ± • 

Proof. Since the proof of this proposition only depends on the structure of the kinetic equation instead of 
the particular form of the boundary condition, the details are the same as in Proposition 3.8 in LLS14 . We 
explain the main idea here. The key structure we utilize here is that the coefficients in the added damping 
terms only depends on the average of the damped solution against or ^L~^{pXq)- Hence to remove 

the damping effect, we only need to choose g carefully such that f — g will have zero averages. Bearing this 
in mind, we denote 


t7+= ((mX+, 1 ,/), ■ 

Uo = 5 • • 

Uc,0=(^(^h‘^ ^(M^0,i)i/^ ; (('^’1+^^)^ ^{{vi +u)XoAXQ^,ygA)y'j 


U- = [{gX_^uf), ■■■, 

h'^0,uo ’ ’ 


and 


Uf = 




ul, uA A 


C,0 


By multiplying X+j,X-^i,Xo^kA ^{xiXo,m) to (3.1) and integrating over v G K“, we have 


do,U + AU = 0, 


(3.17) 


(3.18) 


where the coefficient matrix A is 


A = 


where D± are positive diagonal matrices and 
Bij = (^gXg^i, C ^(/iAoj)^ 


aD+ 

—aD- 

0 

0:^21 

0:^22 

0 

0 

aB 



I + aB 

aD 


(3.19) 


7+X70 


^22,jk — B ^(^Ao_fc)^^ 


7- X7o 


Dij = (gC ^(^Ao,i), C ^{gXoj) 


70x70 \ / 70x70 

where B is symmetric positive definite and D is symmetric. Thus A is a matrix of size ( 7 + + 7 _ + 270 ) x 
( 74 - + 7 _ + 270 ). The proof of Proposition 3.8 in LLS14 shows that A has 7 _ + 70 negative eigenvalues 
Moreover, A is of rank 7 + + 70 since the original kinetic equation (1.2) satisfies the uniqueness 
property in Lemma 2.2 Note that by the boundedness, all the solutions to the damped equation (3.1) will 
be orthogonal to span{vi}7rj^^°. Hence for any solution / to the damped equation, there exists a unique set 
of {co,i}7=i U {c+,j}l=i 9 defined in (3.15) with these coefficients, we have 

AtJ, = AUf , 


which is equivalent to (3.16). 


□ 


Now we can construct the solution to the original kinetic equation. 
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Proposition 3.5 Let ’s and ’s be the coefficient of g given in Proposition 

vq 1^+ 

^ = / - 5 + X! + X! ^+ 3^+3 ■ 


3.5 Let 


Proposition 3.6. Let f be the solutions to the damped equation (3.11 with h and g the function defined in 

(3.20) 


i=l 


i=l 


Then fj is the unique solution to the original half-space equation (1.2). 
Proof. Note that pi = f — g satisfies 
p-dxffi + = 0 : 


Pi\ 


Ai >0 


/ 1^0 


^ I /i .>0 


\i=l 


i-1 


771 0 , 


a; = 0 , /I > 0 , 


as a; —)■ 00 . 


Hence fj defined in (3.20) satisfies 

pLdxfj + Cfj=Q, 
fj = h + ICfj, 


>'+ 


P 


*^+3^+3 ’ 


a: = 0,71 > 0 , 


as X —>■ c». 


i=l 


i=l 


where co,i,co,j are the coefficients defined in Proposition 3.5 We thereby have recovered fj as the unique 
solution to ( 1 . 2 ). □ 

Combining the error estimate in Proposition [3]^ and the damping terms, we derive the final error estimate 
for our method as follows: 


Proposition 3.7. Suppose fj is constructed as in (3.20) in Proposition 3.6 with f, g+,i, goj being numerical 
approximations obtained in Proposition 3.3 to the damped equation with appropriate boundary conditions. 

hen the'. 

f - w 


Suppose fh is the unique solution to the equation (1.2). Then there exists a constant Cq such that 

fh-w 


Wfh - p\\r < Co ( ^inf 

^wGTnk 


where ||* ||p is the norm defined in ( |3.2[ ) and 

5nk 


+ inf 
r w^Tnk 


+ ^NK 


(L2(a dn dx))'” 


\m,j - w||r ■ 


wCiTnk ujGFivK 

Proof. The proof of this proposition only depends on the recovery procedures and the quasi-optimality shown 
in Proposition |3.3[ In particular, it does not depend on the specific form of the boundary conditions. Hence 
it is identical to the proof for Proposition 3.9 in LLS14 and we omit the details. □ 

4. Numerical examples 

In this section we show the numerical results of our algorithm for three models, which cover the cases for 
multi-species, multi-dimensional (in the velocity variable), and multi-frequency systems. The three examples 
are: a linear transport equation with two species, linearized BGK/Boltzmann equations with velocity in K^, 
and a linearized transport equation with multi-frequency. We treat these three cases in order. Recall the 
general form of the half-space equation: 

hdxf + Cf = 0, in(0,oo)xV, 


C>« = 3p) + >c{fX^„) 


at X = 0 . 


(4.1) 
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As mentioned before, the boundary conditions for all these examples are either the Dirichlet condition with 
given incoming data or the classical Maxwell boundary condition such that 

/C = Cly/Cy = CXdlQd -\- CXs^S : 

where /Cd is the diffuse reflection and /C^ the specular reflection. For the convenience of numerical computa¬ 
tion, we list out two properties of such /C which can be verified by direct calculation: 


Lemma 4.1. Let K, the Maxwell boundary operator and let /C, K-s, ICd be the operators defined in Lemma 
Then 

(a) = I; 


2.3 


(b) (I + /C) ^ has the explicit form as 


(l + IC)-^ = 


1 


1 H” 


(I - idICd) 


(4.2) 


where jd = C(d{^ + ad + as) ^ 

4.1. Linear Transport Equation with Two Species. 

4.1.1. Formulation . The first example that we consider is the steady radiative transfer equation (RTE) with 


Thomson (Rayleigh) scattering and polarization effect in planar geometry (see Pom73 Section 4.5]). In this 
model, the variables / and Q denote the total intensity and the intensity difference of light. The system 
[Pom73 Eq. (4.211), page 135] depends on the frequency which only serves as a parameter. Hence we simply 


ignore the frequency dependence here. In this case the scattering coefficients ct, CTs in Pom73 are both 


constants. We consider a pure scattering case with no source such that a = Os and rescale a to be one. The 
speed of light c is also normalized to be one. Then the RTE has the form 


p,dxl +1 - 


1 


1 -f 


^P2(ii)P2(m')^ lih-') V - ^ J P2(m) (1 - P2(ll')) Q{t') = 0 , 


pdxQ + <3 — 2 f~2 


y ^(l-p2(At))P2(M')-f(M)d/r' +^ J {^-P2{ij))[l-P2{F))Q{T)<^T^ =0, 

where P 2 {p) = ~ i i® second-order Legendre polynomial. 

4.1.2. Properties of C and 1C. Denote / = (7,(3)^. Then the collision operator C has the form 

{1+ Ip2{p)p2{p')) -i(l-P2(M))P2(/i') 


(4.3) 


(/:/)(/r)=/(/r)-(E(Ai, •)/(•) 


S(li,/i')=ri 

\-UI-P2{t'))P2{p) \{1-P2{p)){l-P2{p')), 


(4.4) 

where we recall the notation {gi, g 2 )^ gig 2 d/i. In this case, we have dcr = d/i. First we check that 

Lemma 4.2. The scattering operator C defined in (|4.4)) satisfies (PLIJ-(PL5J. 


Proof. The attenuation coefficient a in this case is a{pi) = 1. One can then directly check that C is self-adjoint 
and Null£ = span{Ao} with 

Xo = (1,0)^ and C-\pXo) = g-Xo . (4.5) 

Property (PL2) is also readily verified by the boundedness of E(/i,/i'). Furthermore, we can show by direct 
calculation again that there exists a constant /32 > 0 such that 


f ■ Cfdg> P 2 


V^f 




(4.6) 


where is the projection onto (Null£)■*'. Hence (PL4) holds. 
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To show that C satisfies (PL5), we prove a more general statement: suppose C satisfies (PL1)-(PL3) 
together with (4.6 1 , then there exist constants a,ao >0 such that (PL5) hoids. Indeed, write f = + 

fo + f++ f-, where 




/o = ^(Xo..,/) Xo,^, /± = E(^±.^>/) - 


Recall that XQ^i,X±j are defined in Section [O] Then by Cauchy-Schwarz there exist ci,C 2 > 0 such that 


f,Cf) >132 


(L2(d<T))" 


^(/xXo,.,/) f) > 0 , 


^(/r£ ^{fiXo,i), f) > 


Cl 


fo 


/) + /) > 


L'^(da) 
2 


C2 


L 2 (do-) 


/+ 


L'^(da) 


Cl 


/+ 


L2(do-) 


/- 


L2(dcr) 


C2 


/- 


2 

L 2 (dCT) 

2 

L2(dcr) 


Therefore, if the coefficient a in the definition of £d in (2.4) is small enough, then 


/• Cdfda > ^ 

z 


(L2(d^))" 

Furthermore, 

[ /•/Id/dcr > a V /) 

Jy ^ ^ 


CiQ; 


/+ 


(L 2 (da))’' 


/- 


(L2(da))" 


(4.7) 


> cio; 


/o 


(L2(d<T))" 


C 2 a 




(L2(d<T))" 


/+ 


(L 2 (d^))" 


/- 


(i"(dcT))" 


(4.8) 


Hence by multiplying (4.7) by max{^, we have 


/ • £d/ dcr > (Ti 


(L2(da))- 

for some cti > 0 which depends on a. Applying the above estimates to £ given in ( |4.4[ ) (note that in this 
case we only have 77° or Xq), we conclude that such £ satisfies assumptions (PL1)-(PL5). □ 


The boundary operator K is defined as 

«(/£<.) =”A£+»a.(/1 „J 


(4.9) 


= 


V '<0 




Xn 


J^,^^H'Xo-Xodf,' 




where Xg is given in (4.5), and ad, cis > 0, 0 < cid + < 1. We have 


Lemma 4.3. The boundary opeartor X defined in (4.9) satisfies (fPKJ. 

Proof. Since Lemma |2.1| only covers the scalar case, we show the details of proof for the current vector 


case. As commented in Remark 2.1 we only need to show that Xd satisfies (PK). To this end, we use the 
symmetry of Xq in p and Cauchy-Schwarz to obtain that 


p\Xdf\‘^dp = 


' fi>0 


f f 



(/ ^ 

^0 

d/r) 

\d/ 4>0 


/ 


Ii.<o\T\f-^odp 


// 4 >o ■ Xo dp 


< 


Mfl'^dp, 


’ Ai <0 


which is the desired property (PK). 


□ 
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The above two Lemmas show that the theory in Section [^applies to equation (4.31. The details for terms 
in the weak formulation are as follows. First, 


^df = Cf - aiiXo fj- a^J.C (pXo) {jiC (/rXo), / 


= Cf- / 


The Maxwell boundary condition is (recall (4.9)) 




(4.10) 


Thus, by Lemma 4.1 the boundary operators in the weak formulation /) = l{(j)) are 


(1 + ^)-' = 1C,) , 

1 + <Ts \ 1 + Cld + Q!s J 

‘ 2 ad 


(1 + /C)-^(I-/C) = 


1 - Gs 


1 - 


1 + Os \ (1 ~ Cls)(l + Cld + Cls) 


ICd 


where 


lCdf = ICd{I, Qf = 2([ ^7(0, /i') d^') (1,0)^ . 

\Jd .’>0 / 

Let go be the special solution to the damped equation with the boundary condition 

5oUo = (^0 -/C(Xo|^^q)) +/C(go|^<o)’ ^ = ^- 


(4.11) 


Then the true solution is given by 


fh = f — Chgo + ChXo . 


where / is the solution to the damped equation with the boundary condition (4.10) and c/j = 


{d-Xo, go ) 


4.1.3. Algorithm. For this example, we choose half-space Lagendre polynomials as the basis functions for 
each component of Namely, we first find the half-space Lagendre polynomials by: 

/ ^rn{fj‘)^n{g) dg. = 6mn , m,n = l,2, •••. 

JO 

Then we use the even-odd extension to obtain the basis functions for the finite-dimensional space over [—1,1]: 

-pi r J, I 2A^— 1 

Tat,! = span{ 0 ^}^^^ , 


where 




4 'mi.g') 1 /7 s [ 0 , 1 ] , 

(fmi-g), /re [- 1 , 0 ) 


and (j) 2 m-i{g) = 


4^m{g )) M G [ 0 , 1 ], 

/rG[-l,0). 


In this multi-species case, / is a two dimensional vector and the basis function for / is chosen to be: 

1 4N-2 

Fjv = span ■ 


>1 4iV —2 

4>m f 

J m=l 


with 

*“=( 0 

Using these basis functions, the ODE system becomes 


and (j)m+ 2 N-i = 


, m = l,... ,47-2. 

(Pm{x) I 


A—d = Bd, 
da; 
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Figure 1. (Example 4.1) For the two-species RTE (4.31, we set h = (2fj., and ad = 
Us = 0. The left panel shows the numerical solutions at a; = 0 for both species using 31 (i.e. 
N = 16) basis functions. The right panel shows the convergence rate of the Q component 
of the RTE. 


where 


^mn — / /10m ‘ 0n d/l , 


Bmn — / 0m ‘ d/l 


J-1 J-1 

for m, n = 1, • • • , AN — 2. Therefore both A and B are of size {AN — 2) x (4A^ — 2). 

4.1.4. Numerical Results . In this part we show the numerical results regarding this multi-species model with 
pure incoming data. 

Example 4.1 (Incoming boundary condition). Set h = (2/r,/i)^ and ad = as = 0. The numerical solutions 
at a; = 0 for both components are shown on the left in Figure The plot on the right in Figure shows 
the convergence rate of the second component Q where we observe an algebraic convergence rate. This is 
within expectation, as even though the even-odd decomposition captures the jump discontinuity at /a. = 0, 


the solution still has a weak derivative discontinuity near /a = 0 CLT14 ,TF13 

4.2. Linearized BGK Equation. 


4.2.1. Formulation . The second example we consider is the time-independent linearized BGK equation for 
a single species with its velocity v = {fj,,Vy) G Here we normalize the wall temperature and denote the 
absolute Maxwellian M as the wall Maxwellian such that 

1 

M{v) = -e-^. 

Suppose F is the density function for the nonlinear BGK equation. Let / be the perturbation such that 

F = M + Mf. 

Then the linearized collision operator C has the form 

Cf = f-Vf, (4.12) 

where V : L‘^{Mdv) —)■ span{l,a;, |up} is the projection operator. In this case, du is the usual Lebesgue 
measure and dcr = M d?;. 
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The boundary operator K. is given by 

ICf = adJCdf + asJCs (/|^<o) 

= a. \AM Md.) + »./U.(RU ,4.13) 

= adV^ [ Im'|/(w) Mdi; + aj\ (JZv ), 

J /i '<0 

where TZv = V 2 ), ad, Os > 0 and ad + as < 1. We have 


Lemma 4.4. The linear operator C defined in (4.121 satisfies (P'L1)-(P'L5) and the boundary opeartor K, 
defined in (4.13) satisfies (T*KJ. 


Proof. It is classical to show that C satisfies (PL1)-(PL4) and also the stronger coerciveness (4.6). Therefore 
by the proof of Lemma 4.2 condition (PL5) also holds. The boundary operator ICd in this case fits the form 
in Lemma 2.1 which guarantees that (PK) holds. □ 


Remark 4.1. Our method can also be applied to the linearized BGK equation for multi-species. The lin¬ 
earization of these species is chosen slightly differently depending on whether there is diffusion reflection or 
not. In the case where there is nontrivial diffuse reflection from the wall, the equilibrium state for different 
particles will all be the same, which is the Maxwellian M given by the wall. We can then linearize the 
vector-valued density function F = (Fi, • • • , Fm) as 

F = M(1,... ,l)^ + M(/i,.-- ,Uf. 


On the other hand, if there is only the incoming data and/or the specular reflection, i.e., ad = 0, then we 
allow the equilibrium states of different species to be different. In this case, we linearize F = {Fi, ■ ■ ■ ,Fm) 
as 


Fi = Mi + \fWifi , f = 1, • • • , m . 

The advantage of this linearization is that the function space for / is given by (L^(di;))'" instead of the 
weighted-L^ by various Maxwellians Mi for each component fi. 


To define the particular damped operator in the case of linearized BGK equation, we compute the eigen- 
modes in Null£: 

= /^^ + ■Py ~ 4 , Xo,2 = Vy , 

= fjf Vy 2p ,, X— = pfi + Uy — 2/i. 

The associated eigenspaces are 

= span{Xo,i, Xo, 2 } , = span{X±}. 

Moreover, is computed as 

C ^(/iXop) = pXop = -f Ily — 4), C ^{pXop) = pXo,2 = ■ 

Hence, the damped operator has the form 

Ldf =Cf + apX+ {pX+, /)„ + apX_ {pX_, f)^ 


+ a'£,l‘^uAXu.f),+«Y, 






where (/i, /2 = /]^2 fif 2 ^ du. The boundary condition is given as 

/l4>0 = '' + '=(/U) . 


X = 0 , 
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where recall that JC is given in (4.131. By Lemma 4.1 the boundary operators in the weak formulation are 

1 


(i+/cr^ = - 
(i+^)-'(i-^) = 


Cts 


I- 


ad 

l^d I , 


1 + CKd + Of 

2 ad 


(1 - as)(l + ad + as) 


1C, 


where Kdf = /^,>o ^z;. 


In order to obtain the original solution to (1.2), we construct the special solutions 30 , 1 ) 50 , 2 , 9 + such that 
they satisfy the damped equation and the boundary conditions respectively: 

50,1 |^>o = (^0,1 -/C(^o,i|^<o)) +^(5o,i|^<o), 

50,2|^>0 = (^0,2 -/C(Xo,2|^<o)) +^(50,2 |^<o) , 

5+Im>o = (^+ “^(^+L<o)) +^(5+|^<o) ■ 

The true solution is then given as 

fh = f — (co,l 50 ,l + Co, 290 ,2 + C+ 9 +) + (co,i-^o,i + Co,2-^0,2 + C+X+) , 
where the coefficients co,i,co,2,c+ satisfy that 

(M-^O,!) g) = (5-^0,1, /) ) (5-^0,2) 9) = (5^0,2, /) , (5-’^+) 9) = (5-^+) /) 


with g = co,i 5 o,i + co, 290,2 + c+ 5 +. The unique solvability of 9 is guaranteed by Proposition 3.5 


4.2.2. Algorithm: For the 2D-BGK case, we build the basis functions upon Hermite polynomials. Since the 
solution is regular in Vy, we use the full Hermite polynomials on M for Vy. To take into account of the 
jump discontinuity in 3 = Vx, we apply even/odd extensions of half-space Hermite polynomials on [0,oo). 
Specifically, the half-space Hermite polynomials satisfy 

/ Bjji{^fd)I3ji(^ld)€ ^ d^ — ^mn • 


Performing the even and odd extensions of gives 

H„(^)/v^, 5>0) 




Bn{-g)/'/2, 3<0 


and B^{p) = 


Bn{td)/V2, 

- Bn{-n) / V 2 , 


5 > 0 , 

fi < 0 . 


Then the set of basis functions for the finite dimensional space in g is given by 


r.,A = {<^., 2 „-i};L"'i' U {<^., 2 „Ci = U 


O 1^+1, 


> N 


=1 


The set of basis functions in Vy is 


r 2 /,Ar = ,where 




The basis for the approximation solution /tv is then expanded by r 2 ;,Ar 0 ^y,N such that 

2iV+l N 

fNih,Vy) = EE 


m—1 n—1 


The ODE system still has the form 


A—a = Ba , 
ax 
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d-/i I Spq ^ ^mn — / (^a:,n^y,g) ] (m? ) d-//d'Uy . 


where 

Apq — 

n 


4.2.3. Numerical Results . Examples are shown in Figure and Figure for both the pure incoming data 
and the Maxwell cases. 

Example 4.2.1 Incoming boundary condition. In the first example, we use incoming boundary condition. 
Figure [^verihes that if h S 0 then the solution is simply = h, as expected from the theory. 

Example 4.2.2 Maxwell boundary condition. In the second example show in Figure]^ we set the accommo¬ 
dation coefficients to be (02,0:3) = (0.3,0.4). This time if h is chosen such that 


h = (A: = 1,2) or 


h = X: - IC{X. 


I y<0 


then the solution is fh = ^o,fc or which again is consistent with the theory. 

4.3. Multi-frequency linearized transport equation. 

4.3.1. Formulation . In this third example, we consider a linearized BGK-type of equation that models 
phonons with a continuous range of frequencies MCMYll . We consider the time-independent half-space 
equation. Let w G (0,ojm) be the angular frequency of phonons and F be the density function such that 
F = F(x,fi,uj). 

The stationary nonlinear equation has the form 

— Fbe 


^v{uj)dxF = -- 


t{uj) 


(4.14) 


where v{ui) is the group velocity and t{uj) is the relaxation time. They are both frequency-dependent. We 
assume that v{ui) > 0, although it may not have a positive lower bound. The equilibrium state Fbe is given 
as the Bose-Einstein distribution function such that 

1 


FBE{i.^, T) = 


(4.15) 


where h is the reduced Planck constant, kB is the Boltzmann constant, and T is the temperature. Given a 
reference temperature Tg, we linearized F{u!,T) around F{uj,To) such that 

F(w,r) = F(w,ro) + G(w,T). (4.16) 

The resulting equation has the form 

G-C^AT 


pLv{(jj)dxG = -- 


r(w) 


(4.17) 


nw/(fc^To) 


where C^i = ^^h^/{kBTo)_iy mass conservation, we have 


r-" r G 

lo J-1 t{uj) 


d/r du! = 


which gives 


AT= ^ 


G 


©0 Jo J-I 'r(w) 


d/i dw 


r- r 
'0 J-l t{uj) 


with 0n = 


d/i dw ) AT, 


a 


1-1 t(w) 


d/i dw. 


Equation (4.17) then has the form 


/iu(w)5,G = -^ P r -^d/idc^) . 

t{u}) \ ©0 Jo J-I t{uj) J 


(4.18) 
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Figure 2. (Example 4.2.1) These four rows of figures demonstrate the results computed 
using h = Xq i, h = Xq 2 j h = X^ and h = X_, and the three columns show h, recovered 
solution fh and the difference h — fh respectively. For the first three cases, the solutions 
satisfy that fh = h since h G 0 i? +. Solution to the last case does not satisfy fh = h 

since h G H~. In all examples, we use 31 basis functions along each direction. 
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Figure 3. (Example 4.2.2) These four rows show numerical results using boundary con¬ 
ditions provided hy h = X — where X = X-,Xo,i,Xo, 2 ,^+ respectively. The 

accommodation coefficients are set as {ad,cts) = (0.3,0.4). The four columns (from left to 
right) are: h, X, recovered result ft, and the difference ft — X. In the bottom three cases, 
one recovers Xo,i,Xo, 2 ,X+ as the solutions as expected. In all examples, we use 31 basis 
functions along each direction. 


The operator on the right-hand side of (4.18) is not self-adjoint in the space L^(d^dw). However, we show 
below that it is symmetrizable. Indeed, if we define 

a 


Pu, = v{uj)t{lu) , 


G = 




/, 


(4.19) 
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then / satisfies 


= -^f + 


Pui Jo J-1 yfRj 


^ d. ^ ^ 




= -m): 


(4.20) 


where dcr = d^dw is a probability measure by the definition of ©o- Note that the linear operator C 

on the right-hand side of (4.20) is now self-adjoint. 

In the previous two examples, both linear operators C satisfy the classical “coercivity” condition (4.6). 
However, this property ceases to hold in the current example. Nevertheless, we show that with the help of 
the added damping terms, condition (PL5) is still true. 

Lemma 4.5. Suppose 0 < w < uim, 0 < Pui < bo, and f = f{pL,ui). Let L he the scattering operator given by 

rl f 




Pbj Jo J-i y/J^ 


da, 


(4.21) 


where dtr = d^dw is a probability measure in (/r,a;). Then 

(a) C is self-adjoint, nonnegative, and Null£ = span{v(S(0}. Moreover, y/]5fj G ■ 

(b) Denote Xq = jfXd __ ^ D^ji^ie 

y/JvP<^ d(T 

Cdf = Cf + a(^fiXo{pXo, f)^^^ + pC-\^iXo){^iC-\yiXo), . 


(4.22) 


where {gi, fli 9 iff 2 da. Then for a > 0 small enough, there exists a^ (depending on 

a) such that 


j^fCdfda > (72 ||/|li2(d„) . 
Hence C satisfies all the assumptions (PL1)-(PL5). 


(4.23) 


Proof. By the definition of £, we have 

rOJm P 1 


fCfda = 


/o 


LUXiwX 

f' f f / / \ " 


dcr 


a/ Pui \ Pui j , 


dcr > 0 . 


(4.24) 


This shows £ is nonnegative and Null£ = spanjy^}. By direct calculation we have = 0. Hence 

culation that 

(m(v^)^) 


G Denoting Xq — = , one can verify by direct calculation that 

viv 


L-\pX^) = 


fy/3^da 


Hence, 


p.£ (pXo), f 




l/iyAhio- 

1 

/v dcT 


(fff / / / \ 


y/1^ \ yfWJ / 


ao 




yfKj 


fl^UJ 
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where ao = > 0. Thus, 

V /v /3„ do- 


/O d-l 


/£,/dcT = ifCf)^^^ + a ifiXo, + a (m£-'(m^o), /)',^ 


2 / r \ ^ 

aao / / \ 


a 


,2o2 / 






/ 


/ 


/-iIva: \yA:/.,J 2°\vs; 


/i,u; / \ V / //,a; 

where the last inequality follows from choosing a small enough and then applying Cauchy-Schwartz and 

1^—I 2 

(4.24). Let a = min{ Then by the assumption that (3^ G (0, 6o]) we have 


fCdfda > a 


/o J-i 


f 




> Cr2||/|li2(do) > '^2 = —, 

L2(do) ^0 


which proves the coercivity of the damped operator Cd on (dtr). Hence C satisfies all the assumptions 
(PL1)-(PL5). □ 


The type of boundary conditions we use here is the case where the wall does not change the frequency 
of the phonon. More precisely, for each u; > 0, the boundary condition for the original nonlinear equation 
reads 

Fbe(w, To) 

IJ,\v\ujjr -j - — 

f fi<0 


^^\n>0=^d [ \^1\V{U})F d^lj- 

J tl<.0 J n 


V>0 l/^k(w)TBE(w,To)dAi 

= 2adJ \^^\v{uJ)Fd^J. + asF\^^^(-n,u). 




Linearizing F as in (4.16) and (4.19), we obtain the linearized boundary operator K, as 


Kf = adfCdf + asICsifl^^g) = 2ad J |li|/d/r -f Q^s/|^^g(7^^)), 

where ad, as > 0, and + Os < 1. Now we verify that 

Lemma 4.6. The boundary operator K, defined in ( |4.25| ) satisfies (T’KJ. 


(4.25) 


Proof. Again we only need to show that ICd satisfies (PK). By the definition of ICd in (4.25), we have 

[ Ai|/Cd/p dcr = 4 / /if"/ \p.'\f{p.',uj)dp'\ k do; dp, 

Jfi>o \Jfi'<o J Co'P(w) 


< 2 






a 


[ lAi|/^(M,w)dCT, 

' /i<0 


0ot(w) 


du)dp! I dpL 


which shows (PK) holds. 


□ 


In summary, the damped equation has the form 
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The boundary condition is given as 


f\^^o = h + JCf, 


(4.27) 


where K is given in (4.251. We again have 


(I + /C)-^ = -— 

1 + a, 

1 + Gs V 


OLd -p 


1 + Gd + Qfs 

“^ad 


where 


(1 - as)(l + ad + as) 
Kdf = 2 1 /r/d/r. 


1C, 


(4.28) 


^ M>0 

The special solution go is constructed as 

fJ-da:go + Cdgo = 0 , 

5o|^>o = (-^o-^(^o|^<o)) +^(5o|^<o)’ M>0. 

Finally, the exact solution for equation ( |4.20| ) with boundary condition (4.25) is 

fh = f — Chgo + ChXo , 

where / solves the damped equation (|4.26|) with the boundary condition (|4.27|) and Ch = ^ 

(/rAo, go)^^sj 

4.3.2. Numerical Results. We discretize the w-variable uniformly and replace the integral in w by the trape¬ 
zoidal rule. The resulting system can be viewed as a multi-species system. Hence the construction of basis 
functions is the same as in Section 14.1.31 

We again show examples with both pure incoming data and Maxwell boundary condition. For computa¬ 
tional convenience, we modify the equilibrium state as 

Fb£;(w,T) 

Example 4.3.1 Incoming boundary condition. In the first example, we set 

C I 

—rp = w exp (— w/lOOO) , t{ijj)v{uj) = — , a;G[I,8]. 

T\Uj) UJ 

Then Null£ = spanji/T^ w)u(a;)} = spanly'T/o;}. Figure shows that if /i = Xq, then the numerical 
solution is in good agreement with the analytical solution where fh = Xq. 

Example 4.3.2 Maxwell boundary condition. In the second example, we take the same Cu,, t{uj) and v{uj) 
as in the previous example and set the accommodation coefficients as (ad, as) = (0.3,0.4). Once again if 
h = Xq — /C(Xo|^^p), the exact solution must be fh = Xq. The numerical solution demonstrated in Figurej^ 
shows a good match with the exact solution. 
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Figure 4. (Example 4.3.1) The two figures on the left are boundary data and the difference 
{fh{x = 0) — h) obtained for h = \/l7w S Null£. The two figures on the right are the 
solution at the boundary and the difference (/^(O,/i, w) — h) obtained for h = In this 
case we expect the difference to be of order 0{1). In both cases, we use 31 basis functions 
in ^ direction and sample 8 grid points along lo. 



Figure 5. (Example 4.3.2) In the second row we use h = Xq — as the incoming 

Dirichlet data. The four plots show h, Xq, numerical result fh and the recovery difference 
fh — Xq- The accommodation coefficients are set as {ad,as) = (0.3,0.4). We recover the 
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